# This code replicates Figure 3a
rm(list=ls())
setwd("")
library(multiwayvcov)
library(foreign)
library(ggplot2)
library(lme4)
data <- read.dta("data.dta")

# Write function to extract coeficients
get_coef <- function(sm, vars = 1:4) {
  coefs <- sm[vars, 1]
  se <- sm[vars, 2]
  names <- row.names(sm)[vars]
  data.frame(var = names, coef = coefs, se = se)
}

# Write function to estimate cluster-robust standard errors (robustness)
cluster <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm), 2 , function(x) tapply(x, cluster, sum, na.rm=T));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) 
}

# Function to estimate models
get_robust_coef <- function(dv, ivs, df, fe = T) {
  
  # get formula
  formula <- paste0(dv, " ~ ", paste(ivs, collapse = "+"))
  
  # fe adjustment
  if (fe == T) formula <- as.formula(paste0(formula, " + factor(village) - 1"))
  if (fe == F) formula <- as.formula(formula)
  
  ## get df, do model
  df2 <- df[complete.cases(df[, c(ivs, dv)]), ]
  if (fe) {
    model <- summary(lm(formula, data = df2))$coefficients
  } else {
    model <- cluster(df2, lm(formula, data = df2), df2$village)
  }
  
  ## get index of coefs to pull out
  v <- if (fe) 1:4 else 2:5
  c <- get_coef(model, vars = v)
  c$outcome <- dv
  c$method <- ifelse(fe, "Fixed Effects", "Robust SE")
  c
}

# Define treatment variables and control variables 
ivs <- c("treatment_adventist", "treatment_mixed", "treatment_peruana", "treatment_maranatha")
ivs_ctrl <- c("treatment_adventist", "treatment_mixed", "treatment_peruana", "treatment_maranatha", "age", "male")

df <- data

# Estimate models, extract coeficients
res_obed <- lapply(c("exp1_leave", "exp1_leave_state", "exp1_leave_church"), get_robust_coef, fe = F, ivs = ivs, df = data)
res_obed_ctrl <- lapply(c("exp1_leave", "exp1_leave_state", "exp1_leave_church"), get_robust_coef, fe = F, ivs = ivs_ctrl, df = data)

# Combine coefficients
res1 <- rbind(do.call(rbind, res_obed), do.call(rbind, res_obed_ctrl))

# Estimate linear models to double check and to then estimate wild cluster bootstrap-T standard errors
m1 <- lm(exp1_leave ~ treatment_adventist + treatment_mixed + treatment_peruana + treatment_maranatha, data=data)
m2 <- lm(exp1_leave_state ~ treatment_adventist + treatment_mixed + treatment_peruana + treatment_maranatha, data=data)
m3 <- lm(exp1_leave_church ~ treatment_adventist + treatment_mixed + treatment_peruana + treatment_maranatha, data=data)
m4 <- lm(exp1_leave ~ treatment_adventist + treatment_mixed + treatment_peruana + treatment_maranatha + age + male, data=data)
m5 <- lm(exp1_leave_state ~ treatment_adventist + treatment_mixed + treatment_peruana + treatment_maranatha + age + male, data=data)
m6 <- lm(exp1_leave_church ~ treatment_adventist + treatment_mixed + treatment_peruana + treatment_maranatha + age + male, data=data)

# Estimate wild cluster bootstrap-T standard errors
boot_firm1 <- cluster.boot(m1, data$village, boot_type = "wild",wild_type = function() sample(c(-1, 1), 1))
boot_firm2 <- cluster.boot(m2, data$village, boot_type = "wild",wild_type = function() sample(c(-1, 1), 1))
boot_firm3 <- cluster.boot(m3, data$village, boot_type = "wild",wild_type = function() sample(c(-1, 1), 1))
boot_firm4 <- cluster.boot(m4, data$village, boot_type = "wild",wild_type = function() sample(c(-1, 1), 1))
boot_firm5 <- cluster.boot(m5, data$village, boot_type = "wild",wild_type = function() sample(c(-1, 1), 1))
boot_firm6 <- cluster.boot(m6, data$village, boot_type = "wild",wild_type = function() sample(c(-1, 1), 1))

# extract SEs
wildse <- rbind(coeftest(m1, boot_firm1)[7], coeftest(m1, boot_firm1)[8], coeftest(m1, boot_firm1)[9], coeftest(m1, boot_firm1)[10],
                coeftest(m2, boot_firm2)[7], coeftest(m2, boot_firm2)[8], coeftest(m2, boot_firm2)[9], coeftest(m2, boot_firm2)[10],
                coeftest(m3, boot_firm3)[7], coeftest(m3, boot_firm3)[8], coeftest(m3, boot_firm3)[9], coeftest(m3, boot_firm3)[10])
res_obed <- cbind(do.call(rbind, res_obed),wildse)
wildse <- rbind(coeftest(m4, boot_firm4)[9], coeftest(m4, boot_firm4)[10], coeftest(m4, boot_firm4)[11], coeftest(m4, boot_firm4)[12],
                coeftest(m5, boot_firm5)[9], coeftest(m5, boot_firm5)[10], coeftest(m5, boot_firm5)[11], coeftest(m5, boot_firm5)[12],
                coeftest(m6, boot_firm6)[9], coeftest(m6, boot_firm6)[10], coeftest(m6, boot_firm6)[11], coeftest(m6, boot_firm6)[12])
res_obed_ctrl <- cbind(do.call(rbind, res_obed_ctrl),wildse)
res1 <- rbind(res_obed, res_obed_ctrl)

## Creating coefficient plot

# Renaming
res1$var <- ordered(res1$var, levels = levels(res1$var)[c(2, 4, 3, 1)])
lvls <- c("Maranatha", "Peruana", "Mixed", "Adventist")
levels(res1$var) <- lvls
res1$outcome <- rep(c("Obedient to authority combined", "Obedient to state authority", "Obedient to church authority"), each = c(4, 4, 4, 4, 4, 4))

# Plot
pd <- position_dodge(0.4)
res1$group <- c(rep("With controls",12), rep("Without controls",12))
p1 <- ggplot(res1, aes(x = var, y = coef)) + geom_hline(yintercept = 0, linetype = "dashed")
p1 <- p1 + geom_errorbar(aes(ymin = coef - 1.96*wildse, ymax = coef + 1.96*wildse, group = group), position = pd, width = 0.1)
p1 <- p1 + geom_point(aes(shape = group), position = pd, fill = "white", size = 3)
p1 <- p1 + coord_flip() + facet_wrap(~ outcome) + theme_bw()
p1 <- p1 + scale_shape_manual(values = c(21, 22)) + ylim(-0.5, 0.5)
p1 <- p1 + xlab("") + ylab("") + theme(legend.position = "bottom", legend.title = element_blank())
p1

# Save plot
ggsave("01_obedience_exp.pdf", width = 10, height = 4)

